Fluid dynamics of bacterial turbulence 



Jorn Dunkel,^ Sebastian Heidenreich,^ Knut Drescher,^ Henricus 
H. Wensink,^ Markus Bar,^ and Raymond E. Goldstein^ 

^DAMTP, Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge CBS OWA, UK 
^ Physikalisch-Technische Bundesanstalt, Abbestr. 2-12, 10587 Berlin, Cermany 
^Departments of Molecular Biology and Mechanical and Aerospace Engineering, 
Princeton University, Princeton, New Jersey 08544 
^ Laboratoire de Physique des Solides, Universite Paris- Sud 11 & CNRS, Bdtiment 510, 91405 Or say Cedex, France 

(Dated: February 22, 2013) 

Self-sustained turbulent structures have been observed in a wide range of living fluids, yet no 
quantitative theory exists to explain their properties. We report experiments on active turbu- 
lence in highly concentrated 3D suspensions of Bacillus subtilis and compare them with a minimal 
fourth-order vector- field theory for incompressible bacterial dynamics. Velocimetry of bacteria and 
surrounding fluid, determined by imaging cells and tracking colloidal tracers, yields consistent re- 
sults for velocity statistics and correlations over two orders of magnitude in kinetic energy, revealing 
a decrease of fluid memory with increasing swimming activity and linear scaling between energy and 
enstrophy. The best-fit model parameters allow for quantitative agreement with experimental data. 



PACS numbers: 87.10.-e,87.10.Ed,87.18.Hf 

A series of experiments over the last decade pTFTO] has 
shed light on generic ordering principles that appear to 
govern collective dynamics of living matter pTHT5] , from 
large-scale animal swarming [Tl |2] to meso-scale turbu- 
lence in microbial suspensions [3]-[8] and micro-scale self- 
organization in motility assays [9l [10]. Although very 
different in size and composition, these systems are of- 
ten jointly termed 'active' fluids, for which there is now 
a range of continuum theories [121 fT4H22] . From these 
have come important qualitative insights into instabil- 
ity mechanisms [TSHUI HOI US] driving dynamical pat- 
tern formation, but a quantitative picture remains in- 
choate; even for the simplest active (e.g., bacterial or 
algal) suspensions uncertainty remains about which hy- 
drodynamic equations and transport coefficients [24l [25] 
provide an adequate minimal description, due in large 
part to the inability of existing data to constrain the 
manifold parameters in these models. One approach to 
remedy this problem is to characterize collective dynam- 
ics as in high Reynolds number fluid turbulence, in terms 
of kinetic energy, enstrophy and spatio-temporal correla- 
tion functions, and to compare with an appropriate long- 
wavelength theory (i.e. Navier- Stokes- type equations). 
We present such an analysis here, measuring collective 
behavior in dense suspensions of the bacterium Bacillus 
subtilis in comparison to predictions of a (fourth-order) 
continuum model for bacterial flow |71 126J . 

Previous experimental studies of bacterial suspensions 
in open droplets d [1 [271 [28] , freestanding films [5l [3 
\25\ [29] , on surfaces (6] [SO] [31], or quasi-2D microfiu- 
idic chambers [7 focused separately on the bacterial and 
fluid components, leaving uncertain how accurately pas- 
sive tracers [32, 33 reflect collective bacterial dynamics. 
The experiments reported here, performed in closed 3D 
microfluidic chambers, allowed near-simultaneous mea- 



surements of cell and tracer motion, and exploit a natu- 
ral reduction in bacterial swimming activity due to oxy- 
gen depletion [8] [271 [31] to obtain data spanning two 
orders of magnitude in fluid kinetic energy. Combined 
with extensive 3D numerical simulations of the model, 
this data allows robust parameter estimates. Quantita- 
tive agreement between experiment and theory suggests 
that this model presents a viable generalization of the 
Navier-Stokes equations to incompressible active fluids. 

Wild-type strain 168 of B. subtilis has cigar-shaped cell 
bodies on average 0.8 /im in diameter and 5 /im long [7]. 
It was streaked on LB medium agar plates from frozen 
stocks. Colonies from these plates were used to inoc- 
ulate overnight cultures in Terrific Broth (TB; Sigma), 
which were back-diluted 1:100 into 100 ml of TB and 
grown to mid-logarithmic phase on an orbital shaker at 
37 °C. These cultures were then concentrated 400 x at 
4,000^ (final volume fraction ^ 50%), and fluorescent 
microspheres (diameter 1 /im, F-8816, Invitrogen) were 
added at a final concentration of ~10^ beads/ml. The 
resulting suspensions were loaded into polydimethylsilox- 
ane (PDMS) microfluidic devices, consisting of a series 
of cylindrical chambers (radius 750 /im, height 80 /im), 
connected by thin channels [71 [36] . The inlet and out- 
let of the device were sealed with vacuum grease, and 
images were acquired in the (x?/)-midplane of the cham- 
bers, ^40 /im above the bottom, using a Zeiss 40 x (NA 
1.3) oil immersion objective and a high-speed camera at 
40fps (Fastcam SA-3, Photron). Movies were recorded 
in pairs for each field of view (768 x 800 pix; 1 pix = 
0.36 X 0.36 /im^), one with bright-field illumination and 
one with fluorescence excitation by a 633 nm laser (B&W 
Tek) at ~20 mW. These movies were taken immediately 
after each other with a ~3min time lag between subse- 
quent pairs. During the ^lOmin imaging period for each 
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FIG. 1: (color online) Flow fields from experiments and sim- 
ulations |35j . (a) Very dense homogeneous suspension of B. 
suhtilis overlaid with the PIV flow field showing collective 
bacterial dynamics. Longest arrows correspond to velocity of 
30 /im/s. (b) Streamlines and normalized vorticity field deter- 
mined from PIV data in (a), (c) Turbulent 'Lagrangian' flow 
of fluorescent tracer particles (false-color) in the same sus- 
pension, obtained by integrating emission signals over 1.5 s. 
(d) Partial snapshot of a 2D slice from a 3D simulation of the 
continuum model (parameters in Table IT]). Scale bars 70 /xm. 



device, the motility of B. suhtilis cells decreased markedly 
due to oxygen depletion [27 . The experimental setup 
yields 2D projected velocities of 3D suspension motion 
(Fig. [T]). Data were analyzed under the assumption that 
the flow structures are isotropic, as verified by test mea- 
surements at different distances from the chamber bot- 
tom. Commercial particle tracking velocimetry (PIV) 
software (Dantec Flow Manager) was used to determine 
the bacterial flow velocity {vx^Vy) from bright-field im- 
ages (Fig. [l^,b), corrected for systematic pixel- locking 
errors Data shown in Figs. |2] and [s] are based on 7 
movie segments (40 fps, each 50 s long) corresponding to 
7 different activity levels. 

Global bacterial flows were quantified by the in-plane 
kinetic energy Exy{t) = {{v^. ^ Vy)/2) and in-plane en- 
strophy Vtz{i) — (^^Z^)' where uj^ = dxVy — dyV^ is the 
vertical component of vorticity and ( • ) is a spatial aver- 
age. While Exy and Vtz fluctuate, their time averages 
{Exy^^z) are approximately constant during the 50 s 
time interval used in the data analysis (Fig. [2]3,c). Over 
two orders of magnitude in energy (Fig. [2]i) we observe 
the linear scaling = E^yjlS?^ with A ^ 24 /im being 
roughly one half of the typical vortex radius. 

Probability distribution functions (PDFs) of the in- 



plane bacterial velocity are approximately Gaussian, 
with a slight broadening due to collective swim- 
ming (Fig. [2|i). The negative values of the equal-time 
spatial velocity correlation function (VCF; Fig. [3|i) indi- 
cate the existence of vortices [7 (Fig. [T]). The VCF is 
remarkably robust with respect to changes in the bac- 
terial activity; in particular, the typical vortex radius 
~ 40 /im, estimated from the first zero of the VCF, 
depends only weakly on the kinetic energy. This result 
is consistent with recent findings by Sokolov and Aran- 
son [8^ for free-standing films. The vortex size in 3D is 
roughly five times larger than for quasi-2D turbulence 
in thin microfluidic chambers [7J, where bacterial swim- 
ming and hydrodynamic interactions are suppressed by 
the nearby no-slip boundaries [36l [37] . Unlike the spa- 
tial VCF, the two-time velocity auto-correlation func- 
tion (VACF) varies systematically with energy or vor- 
ticity (Fig. Isb) , but they collapse when plotted as func- 
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tions of the dimensionless lag-parameter r^J (inset of 
Fig.js})), implying that the higher the activity the shorter 
the memory of the bacterial fluid. Generally, the statis- 
tics of 3D bacterial turbulence differ strongly from con- 
ventional 3D Navier-Stokes turbulence [38ji39j, as bac- 
teria inject energy on the smallest scales, inducing an 
'upward' energy cascade towards larger length scales. 

We infer the flow of the solvent medium from particle 
tracking velocimetry (PTV) analysis of the fluorescence 
images, which only show the tracer particles, assuming 
that they are passively advected. Data shown in Figs. [2] 
and [3] are based on 7 movies (40 fps, length 100 s) at 
different activities. Trajectories of individual tracer par- 
ticles were found with a custom algorithm which, depend- 
ing on seeding density and tracer dynamics, was able to 
identify up to 10^ in-plane tracks, the longest typically 
lasting 5 — 8 s. The effective sample size was insufficient 
to determine reliably the tracer VACFs, but did yield 
global flow properties, velocity histograms and equal- 
time VCFs. The velocity PDFs, calculated directly from 
individual tracer velocities, are approximately Gaussian 
with a peak at small velocities from tracer accumulation 
near the vortex centers (Fig. [2|l). 

Estimates from PTV for the medium VCF and enstro- 
phy were obtained by interpolating tracer velocities on 
a 450 X 450 pix subwindow in the center of the imaging 
plane using MATLAB's Delaunay triangulation with a 
lattice spacing A = 90 Y^pix/A^j, where N j is the mean 
number of tracers detected per frame. The accuracy of 
this reconstruction procedure is controlled by the tracer 
concentration, which was kept low to limit effects on the 
bacteria motion and to avoid tracking ambiguities (typ- 
ically TV/ G [47,144] for data shown in Figs. [2]and[3|. 
As a result, the uncertainties for the PTV data are con- 
siderably larger than for PIV data (see Fig. |2]i). The 
interpolated tracer flow fields were used to estimate the 
kinetic energy E^y^ enstrophy ^z-, and spatial correla- 
tion functions of the in-plane medium flow components. 
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FIG. 2: (color online) Experimental results for bacterial 
and medium fiov^^s, color-coded for activity level, (a) PDFs 
of the Cartesian in-plane velocity components, normalized 
by their mean values and standard deviations, are approxi- 
mately Gaussian (dashed) for both tracers and bacteria, with 
observable systematic deviations. The bacterial flow PDFs 
show slight broadening due to active swimming, which is well- 
reproduced by the model ([T]). By contrast, the PTV distri- 
butions exhibit higher peaks at small velocities due to accu- 
mulation of tracers near vortex centers. (b,c) Mean kinetic 
energy and enstrophy of the in-plane bacterial flow compo- 
nents show moderate temporal fluctuations during the data 
aquisition period, very similar to corresponding PTV data 
(not shown), (d) The time-averaged enstrophy scales linearly 
with the time-averaged energy. Open circles are averages of 
the curves in (b), (c). Errorbars indicate standard deviations. 



vection and nematic interactions, and Ai an active pres- 
sure contribution [26]. For pusher swimmers [36] like 
B. subtilis, general considerations of hydrodynamic [42] 
and nematic stresses [26l [43] suggest that Aq > 1 and 
Ai (Ao — l)/3 > in 3D. The (/3, 'Uo)-terms correspond 
to a quartic Landau-type velocity potential [Ml (TS] [40] 
and are physically motivated by the observation of ex- 
tended jet-like streaming regions in B. subtilis suspen- 
sions at intermediate concentrations [28] . The parameter 
^0 defines the collective speed that would be achieved if 
all bacteria were to move in the same direction. When 
(3^0 the model does not conserve momentum or en- 
ergy, as it describes exclusively the bacterial flow compo- 
nent, which may exchange energy and momentum with 
the solvent. The nonlocal (Fq, r2)-terms encode passive 
and active stresses due to hydrodynamic and steric inter- 
actions. For Ao = 1, Ai = /3 = r2 = and Fq > 0, the 
model reduces to the incompressible Navier-Stokes equa- 
tion. A detailed stability analysis [26^ shows that when 
Ao 7^ 0, /3 > 0, ^0 > 0, F2 > and Fo < this is one 
of the simplest vector models to describe phenomenologi- 
cally the formation of jets and turbulent vortices in quasi- 
incompressible active suspensions. Very recently, the 2D 
version of Eq. ([T]) has been shown to provide a quantita- 
tive mean field description of bacterial meso-scale turbu- 
lence in quasi-2D suspensions [7 . Its applicability to the 
physically more relevant 3D case is first explored here. 

We simulated Eq. ([T]) in 3D with periodic boundary 
conditions using a pseudospectral operator-splitting al- 
gorithm [111 |45] and a pressure correction subroutine to 
ensure incompressibility [3 [26] . Simulation grids ranged 
from 128^ lattice points for parameter pre-screening to 
256^ for statistical analysis. Numerical stability of the 
solver was verified for a wide range of parameters and 



In agreement with the PIV results for the bacterial flow, 
we find again a linear enstrophy-energy relation (Fig.[2]i) 
and comparable vortex radii, using the first zero of VCF 
as an estimate (Fig. |3^). We may therefore conclude 
that, at our very high bacterial concentrations, solvent 
and bacterial flow statistics become tightly linked. 

We now examine how these data compare to predic- 
tions of a theory of active fluids introduced recently 
[3 [26]. This minimal continuum model assumes that, 
at high concentrations, the bacterial flow due to swim- 
ming and advection can be described by a single velocity 
field v{t^ x) and a pressure p(t, x). They obey the incom- 
pressibility condition V • v = and 

ToV\-r2{V^)^v. (1) 

Equation ([l]) extends the incompressible Toner- Tu the- 
ory [TSfHO] with a fourth order term as in the Swift- 
Hohenberg equation [41] . The parameter Ao describes ad- 
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FIG. 3: (color online) Correlation functions for solvent 
(PTV) and bacterial (PIV) flow at different energies and best- 
fit continuum theory (see Table [ij), using the same colors as in 
Fig. 2. (a) Both PIV and PTV data indicate a characteristic 
vortex radius Ry ^ 40 /xm. The decay of spatial correlations 
at small-r depends only weakly on the activity level, (b) Ve- 
locity autocorrelation functions of the bacterial flow collapse 
when the time-lag r is rescaled (inset) by the enstrophy time- 
scale Since oc Exy (Fig. bk), this implies that the 
higher the bacterial activity the shorter the flow memory. 
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FIG. 4: (color online) Isoenergy surfaces (E — 1.7vq) and 
selected stream tubes from the best-fit 3D simulation (visu- 
alized with ParaView) indicate a typical vortex length scale 
and extended band-like regions corresponding to co-aligned 
bacterial jets. See Supplemental Material 35 for a movie. 



space-time discretizations. All simulations were initi- 
ated with randomly chosen velocities. Figure [4] shows 
structure-formation in a typical simulation domain. 

Since in 3D we have Ai (Aq - 1)/3 [26], Eq. ^ has 
essentially five free parameters (Aq, /3, 'Uq, Fq, r2). Two of 
those can be eliminated by choice of appropriate length 
and time units. We adopt a natural unit system such that 
the vortex wave-length scale Ap = 27rY^r2/(— Fq) = 27r 
and vq = 1. In our simulations, the box length is fixed 
as 1/ = 12Ar, corresponding to approximately twice 
the experimental field of view, and the time step as 
At = 0.05Ar/(27r'Uo)- To estimate the three remaining 
parameters (Ao,/3,Fo), we note that Fq and F2 define 
a typical vortex speed Vr = \/— Fq/F2. In the turbu- 
lent regime, it is plausible that Vy is smaller than but 
close to 'Uo, i.e. Vr = C^o where C ^ 1- Furthermore, for 
pushers, the dimensionless parameter Aq should be larger 
than 1, but smaller than for quasi-2D suspensions [7 , 
since nematic (steric) stresses can be more easily avoided 
in 3D; we infer Aq ~ 2. Finally, the acceleration time 
scale To = {/3vq)~^ should be of the order of the vor- 
tex time-scale Ar/Vr- Using these estimates as initial 
values in a systematic parameter scan, and by compar- 
ing with the bacterial PIV data, we obtained the best-fit 
parameters in Table [ij Generally, the VCFs and VACFs 
respond sensitively to parameter variations in the simu- 
lations, suggesting that the estimates in Table |l] are ac- 
curate within 10-15% for quasi- incompressible B. subtilis 
suspensions. As an independent cross-check, we com- 
puted A = {Exy/^zY^'^ from the best-fit simulation us- 
ing Ap ~ 50 /im and found A ~ 29 jiia which compares 
well with the experimental PIV value in Fig. [2]i. We 
stress that the conserved form of the bacterial velocity 
PDFs (Fig. [2^), VCFs, and VAGFs (Fig.|3| implies that 
all our experiments can be fitted by a single set of rescaled 
parameters (Aq, /5, Fq), as it suffices to adjust the physical 



model parameter 


in rescaled units 


in physical units 


Ar - Stti/I 2/(-I 0) 


2tv 


^ 50 /im 




1 


3 — 22 /ym /s 


Ao 


1.7 


1.7 


Vr = y-rg/r2 


0.9 


0.9^;o 


p 


0.1 


1.3 X 10~^(^;oMm)"^ 


Exy 


0.54 


0.54^;g 



TABLE I: Parameters of the best-fit continuum model. To 
match a specific experiment, one must merely adjust the phys- 
ical value of vq by equating Exy — 0.54vo to the corresponding 
kinetic energy value in Fig. |2ji. 



values of vq and Ap to match the kinetic energy and vor- 
tex length at a given bacterial activity level. As evident 
from the flow patterns in Fig. [l] and from the solid curves 
in Figs.|2^ and|3j the best-fit parameters yield good qual- 
itative and quantitative agreement with the experiments. 

For incompressible 'passive' fluids, that are governed 
by the Navier-Stokes equations transport parameters 
have of course been measured for a wide range of ma- 
terials [46 . In contrast, quantitative theories of even the 
simplest active fluids have been lacking. We have shown 
here that the minimal fourth-order vector model [71 [26] 
in Eq. ([T]) reproduces the main statistical features of self- 
sustained 3D bulk turbulence in concentrated bacterial 
suspensions, suggesting that this theory is a viable can- 
didate for the quantitative description of incompressible 
active fluids. Due to the close correlation between bac- 
terial and medium (tracer) flow observed in our experi- 
ments, we expect that this generic model will be useful 
in a wide range of future applications, in particular for 
predicting the effects of confining geometries on collec- 
tive microbial dynamics [471 HH] and for understanding 
the anomalous viscosities of active fluids [24l [25] . 
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